ON SURFACE MESHES INDUCED BY LEVEL SET FUNCTIONS 



MAXIM A. OLSHANSKII*, ARNOLD REUSKENt, AND XIANMIN XU+* 

Abstract. The zero level set of a continuous pieccwisc-affine function with respect to a consistent 
tetrahedral subdivision of a domain in M 3 is a piecewise-planar hyper-surface. We prove that if 
a family of consistent tetrahedral subdivions satisfies the minimum angle condition, then after a 
' simple postprocessing this zero level set becomes a consistent surface triangulation which satisfies 

the maximum angle condition. We treat an application of this result to the numerical solution of 
PDEs posed on surfaces, using a Pi finite clement space on such a surface triangulation. For this finite 
, element space we derive optimal interpolation error bounds. We prove that the diagonally scaled 

mass matrix is well-conditioned, uniformly with respect to h. Furthermore, the issue of conditioning 
of the stiffness matrix is addressed. 
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1. Introduction. Surface triangulations occur in, for example, visualization, 
shape optimization, surface restoration and in applications where differential equa- 
tions posed on surfaces are treated numerically. Hence, properties of surface trian- 
gulations such as shape regularity and angle conditions are of interest. For example, 

^ (-h angle conditions are closely related to approximation properties and stability of cor- 

responding finite elements [1] [2] . 

In this article, we are interested in the properties of a surface triangulation if one 
considers the zero level of a continuous piecewise-affine function with respect to a con- 
sistent tetrahedral subdivision of a domain in M 3 . The zero level of a piecewise-affine 
O^l , function is a picccwisc-planar hyper-surface consisting of triangles and quadrilaterals. 

Each quadrilateral can be divided into two triangles in such a way that the result- 
ing surface triangulation satisfies the following property proved in this paper: if the 
volume tetrahedral subdivision satisfies a minimum angle condition, then the corre- 
£f~) , sponding surface triangulation satisfies a maximum angle condition. We show that 

the maximum angle occuring in the surface triangulation can be bounded by a con- 
, stant max < 7r that depends only on a stability constant for the family of tetrahedral 

■ subdivisions. 

The paper also discusses a few implications of this property for the numerical 
solution of surface partial differential equations. Numerical methods for surface PDEs 
are studied in e.g., [HI SI El El El UH] ■ We derive optimal approximation properties of 
r> ■ Pi finite element functions with respect to the surface triangulation and a uniform 

bound for the condition number of the scaled mass matrix. We also show that the 
condition number of the (scaled) stiffness matrix can be very large and is sensitive 
to the distribution of the vertices of tetrahedra close to the surface. Some numerical 
examples illustrate the analysis of the paper. 

2. Surface meshes induced by regular bulk triangulations. Consider a 
smooth surface T in three dimensional space. For simplicity, we assume that T is 
connected and has no boundary. Let H C R 3 be a bulk domain which contains 
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r. Let {Th}h>o be a family of tetrahedral triangulations of the domain f2. These 
triangulations are assumed to be regular, consistent and stable, cf. [2]. To simplify 
the presentation, we assume that this family of triangulations is quasi-uniform. The 
latter assumption, however, is not essential for our analysis. 

We assume that for each Th an approximation of T, denoted by Th, is given which 
is a connected C 0,1 surface without boundary. In our analysis we assume Th to be 
consistent with Th is the sense as explained in the following definition. 

Definition 2.1. For any tetrahedron S T € % such that m.ea,s 2 (S T H T h ) > 
define T = St nTf,. If every T is a planar, then the surface approximation Th is 
called consistent with the outer triangulation Th ■ 

If Th is consistent with Th, then every segment T = St H Th is either a triangle 
or a quadrilateral. Each quadrilateral segment can be divided into two triangles, so 
we may assume that every T is a triangle. 

Let Th be the set of all triangular segments T, then Th can be decomposed as 



Assumption 2.1. In the remainder of this paper we assume that Th is a connected 
(jO,i sur f ace w ithout boundary that is consistent with the outer triangulation Th- 

The most prominent example of such a surface triangulation is obtained in the 
context of level set techniques. Assume that T is represented as the zero level of a 
level set function <f> and that <ph is a continuous linear finite element approximation 
on the outer tetrahedral triangulation Th- Then if we define Th to be the zero level 
of 4>h then consists of piccewise planar segments and is consistent with Th- As 
an example, consider a sphere T, represented as the zero level of its signed distance 
function. For <j>h we take the piecewise linear nodal interpolation of this distance 
function on a uniform tetrahedral triangulation Th of a domain that contains T . The 
zero level of this interpolant defines Th and is illustrated in Fig. 12.11 



Fig. 2.1. Approximate interface for an example of a sphere, resulting from a coarse tetra- 
hedral triangulation (left) and after one refinement (right). 

In the setting of level set methods, such surface triangulations induced by a finite 
element level set function on a regular outer tetrahedral triangulation are very natural 
and easy to construct. A surface triangulation Th that is consistent with the outer 
triangulation may be the result of another method than the level set method. In the 
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remainder we only need that I\ is consistent with the outer triangulation and not 
that it is generated by a level set technique. 

Note that the triangulation T% is not necessarily regular, i.e. elements T £ 
Th may have very small inner angles and the size of neighboring triangles can vary 
strongly, cf. Fig. 12.11 In the next section we prove that, provided each quadrilateral 
is divided into two triangles properly, the induced surface triangulation is such that 
the maximal angle condition pQ is satisfied. 

3. The maximal angle condition. The family of outer tetrahcdral triangu- 
lations {Th}h>o lS assumed to be regular, i.e., it contains no hanging nodes and the 
following stability property holds: 

sup sup p(S)/r(S) < a<oo, (3.1) 
h>0 seT h 

where p(S) and r(S) are the diameters of the smallest ball that contains S and the 
largest ball contained in S, respectively. The stability property implies that the 
family of tetrahcdral triangulations satisfies a minimum (and thus also maximum) 
angle condition: there exists 6 min > with 

> c(a) > 0, (3.2) 

such that all inner angles of all sides of S G 7ft and all angles between edges of S and 
their opposite side are in the interval [6* m j n ,7r — 8 min ]. The constant c(a) depends only 
on a from (|3.1[) . 

Although the surface mesh T h induced by Th can be highly shape irregular, the 
following lemma shows that a maximum angle property holds. 

Lemma 3.1. Assume an outer triangulation Th from the regular family {Th}h>o 
and let Th be consistent with Th. There exists (f> m in > 0, depending only on a from 
(|3.ip , such that for every S £ Th the following holds: 

a) if T — S (~\Th is a triangular element, then 

< 4>i,T < T - ^min 1 = 1,2,3, (3.3) 

holds, where (/^t ore the inner angles of the element T. 

b) if T = S (~]Th is a quadrilateral element, then 

<t>i,T >0min, i= 1,2,3,4, (3.4) 

holds, where 4>i,T o,re the inner angles of the element T . 
Proof. Let 9 m i n be the minimal angle bound from (|3.2p . Take S e Th- 
We first treat the case where T = S n F^ is a triangle T = BCD, as illustrated 
in Fig. 13.11 Consider the angle := ZBCD. Then either <fi < it — 6 m i n and p.3[) is 
proved with 4> min = 6 m - m or tf> 6 (it — # m in, tt)- Hence, we treat the latter case. Note 
that 

^=sm(ACAF)>sm6 min 
| AC | 

and ZBDC < n-<\> < 9 min < §. Take E on the line through DB such that CE _L DB, 
and F in the plane through ABD such that CF is perpendicular to this plane. Hence, 
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\CF\ <\CE\ holds. Using the sine rule we get 

• / ,,t,™ \AC\ . ,,„ An . \AC\ 1 \CF\ 1 \CE\ 

1 ■sH^BDC)< S ^f^=4^<l. 



sin f min sin tf mi „ sin 



Hence, ZADC < arcsin(^ii^) < 2^%^- holds. This yields 

7 \SlllLin' Sill m in J 



ZADB < ZADC + ZCDB < 2-^— + tt 



With the same arguments we obtain 

ZABD<2 Sin0 
sin 6» min 

Since ZDAB < tt - 6 min and ZDAB = it - {ZADB + ZABD) we get 

mi n< 4-^^ + 2^ -20 =:«?(</»). (3.5) 
sin # min 

Since (f> £ (tt — 8 min , w) C (\tt, tt) it suffices to consider g(<fi) for £ (^tt, tt). Elemen- 
tary computation yields g{\iT) > # m in, g{n) = and g is monotonically decreasing on 
(^7r,7r). Hence the inequality (|3.5|) holds iff < <pQ, where 0o is the unique solution 
in (^7r, tt) of g(4>) = mm . This proves the result in a). 

We now consider the case where T = S PI is a quadrilateral T = ABCD, as 
illustrated in Fig. 13.21 Consider the angle <j> := ZDAB. Then either S (0, # m in) or 
<t> € [# m i n , 7r). We only have to treat the former case. Take E on the line through AB 
such that DE _L AB, and F in the plane through OPQ such that is perpendicular 
to this plane. Hence, \DF\ < \DE\ holds and 

• A \ DE \ 

sin 6 = — — - . 

V AD 
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Furthermore, using ISS- = sin(ZDOF) > sin6* m i n we get 



\OD\ 

DF\ 



S in(ZOAD) = L I ^sm(ZAOD) < j-^J < 

1 \DE\ sin0 
< - = — < 1. 

sin0 min | AD | sin0 min 



This implies 



„ . „ / sin d> \ „ sin < 

ZOAD < arcsin ( ———) < 2— — 



sin « min sin f min 
Hence, since ZDAB = <fi < 2 sin <fi, we obtain 

ZOAB < ZOAD + ZDAB < (l + \ — )2sin0. 

sin 6/ min 

Using ZOAB = tt - ZPAB and ZPAB < ir - ZOPQ < it - 6 min results in 

0mm < (1 + ) 2 Sin <t>- (3-6) 

Sin t^min 

For (j> £ (0, 9 m i n ) the inequality (|3.6[) holds iff > 00: where <fo is the unique solution 

in (0, hir) of 8 min = (l + sin g : )2sin^o- Thus the result in b) holds. 

□ 

The lemma readily yields the following result. 

Theorem 3.2 (maximum angle condition). Consider a regular family of tetrahe- 
dral triangulations {Th}h>o an d a surface triangulation Th = UTeF h T that is consis- 
tent with Th. Assume that any quadrilateral element T = S D Th, S G Th, is divided 
in two triangles by connecting the vertex with largest inner angle with its opposite 
vertex. The resulting surface triangulation satisfies the following maximal angle con- 
dition. There exists </> m i n > depending only on a from (|3.1|) such that: 

< sup 4> t _ T < tt - 0min i = 1, 2, 3, (3.7) 

where ipi.T are the inner angles of the element T . 

Proof. If T = S n T h is a triangle, then ([377]) directly follows from (f3~3"|) . Let 
T = S n T/,, be a quadrilateral, with its four inner angles denoted by 64 > 63 > 82 > 
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6\ > 0. From the result in (|3.4p we have 0j > m i n for all i. The vertex with angle 84 
is connected with the opposite vertex. Let Ti be one of the resulting triangles. One of 
the angles of T\ is 9j with j £ {1, 2, 3}. From 6j > <p m in it follows that the other two 
angles are both bounded by tt — <fi m i n . Furthermore, from 9j = 2tt — 64 — Yl^—x ^ 0i < 
2tt — 9j — 2(f> min it follows that 6j < tt — 4> m i n holds. □ 

In the remainder we assume that quadrilaterals arc subdivided in the way as explained 
in Theorem 13.21 Hence, the inner angles in the surface triangulation Th are bounded 
by a constant 6* < tt that depends only on the stability (close to V) of the outer 
tetrahedral triangulation Th- In particular 6* is independent of h and of how Th 
intersects the outer triangulation Th- 

4. Application in a finite element method. In this section, we use the max- 
imum angle property of the surface triangulation to derive an optimal finite element 
interpolation result. On Th we consider the space of linear finite element functions: 

V h = {v h eC(T h ):v h eV 1 (T) for all T £ T h }. (4.1) 

This finite element space is the same as the one studied by Dziuk in [5], but an 
important difference is that in the approach in [5] the triangulations have to be shape 
regular. In general, the finite element space Vh is different from the surface finite 
element space constructed in [H [9] . 

Below we derive an approximation result for the finite element space Vh- Since 
the discrete surface Th varies with h, we have to explain in which sense is close to 
r. For this we use a standard setting applied in the analysis of discretization methods 
for partial differential equations on surfaces, e.g. [H [5l El [3 [9] . 

Let U := {x £ R 3 | dist(a;, T) < c} be a sufficiently small neighborhood of T. 
We define T h T := {T £ Th \ mcas2(Tn Th) > 0}, i.e., the collection of tetrahedra 
which intersect the discrete surface Th, and assume that T h T C U. Let d be the signed 
distance function to T, with d < in the interior of T, 

d:U^R, \d(x)\ := dist(x,T) for all x £ U. 

Thus r is the zero level set of d. Note that nr = Vd on T. We define n(x) := Vd(x) 
for x £ U. Thus n is the outward pointing normal on T and ||n(x)|| = 1 for all x £ U . 
Here and in the remainder || • || denotes the Euclidean norm on R 3 . We introduce a 
local orthogonal coordinate system by using the projection p : U — > T: 

p(x) = x — d(x)n(x) for all x £ U. 

We assume that the decomposition x = p{x) + d(x)n(x) is unique for all x £ U. Note 
that n(x) = n(p(x)) for all x £ U. For a function v on T, its extension is defined as 

v e {x) := v(p{x)), for all x £ U. (4.2) 

The outward pointing (pieccwisc constant) unit normal on r^, is denoted by n/j. Using 
this local coordinate system we introduce the following assumptions on IV 

p : Th — > T is bijective, (4.3) 

max|d(a;)| < h 2 , (4.4) 

xer h 

max||n(a;) -n h (x)\\ < h, (4.5) 

x£T h 
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where h = sup Te7 -r p(T). In (|4.4[) - (|4.5[) we use the common notation, that the inequal- 
ity holds with a constant independent of h. In (|4.5[) . only x G Th arc considered for 
which iih(x) is well-defined. Using these assumptions, the following result is derived 
in 0. 

Lemma 4.1. For any function u G H 2 {T), we have, {or arbitrary T G J-^ and 
f:=p(T): 

KIIo.t^ Hullo,*, (4.6) 
|u"|i, T ~ |u| 1)f , (4.7) 

Kh.T < M 2 ,T + %ll,T> ( 4 - 8 ) 

where A ~ _B means B < A < B and the constants in the inequalities are independent 
of T and of h . 

4.1. Finite element interpolation error. Based on the results in Lemma |4~T1 
the maximum angle property and the approximation results derived in [I] we easily 
obtain an optimal bound for the interpolation error in the space V^. Consider the 
standard finite element nodal interpolation 1^ : C(r^) — > V^: 

(I h v)(x) = v(x), for all xeV, (4.9) 

with V the set of vertices of the triangles in I\. 
Theorem 4.2. For any u G H 2 (T) we have 

\\u e - I h u e \\ L 2 {rh) < h 2 \\u\\ H 2 {r) , (4.10) 
\\u e -I h u e \\ m{rh) <h\\u\\ H 2 (r) . (4.11) 



Proof. From standard interpolation theory we have 

\\u e - I h U e \\ L 2 (T) </lVkT, 

where the constant in the upper bound is independent of (the shape of) T. Using the 
result in (|4.8[) and summing over T G J- proves the result (|4.10[) . For the interpolation 
error bound in the ii/^-norm we use the results from [T]. For the interpolation error 
bounds derived in that paper the maximum angle property is essential. From [T] we 
get 

\\u e - IhU e \\ H i {T ) < h\\u\\sp(T)- 

Due to the maximum angle property the constant in the upper bound is independent 
of T. Using the results in Lemma 14. II and summing over T G Th we obtain the result 

Amp . □ 

If one considers an i? 1 (r) elliptic partial differential equation on T, the error for its 
finite element discretization in the surface space Vh can be analyzed along the same 
lines as in J5[ . A difference with the planar case is that geometric errors arise due to the 
approximation of r by ITV Using the interpolation error bounds in Theorem 14.21 and 
bounding the geometric errors, with the help of the assumptions (|4.3p - (|4.5|) . results 
in optimal order discretization error bounds. 
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4.2. Conditioning of the mass matrix. Clearly the (strong) shape irregular- 
ity of the surface triangulation will influence the conditioning of the mass and stiffness 
matrices. Let N be the number of vertices in the surface triangulation and {4>i}iLi 
the nodal basis of the finite element space Vt- The mass and stiffness matrices are 
given by 

M = (my)fj =1 , with rriij = / fafads, (4-12) 

A = (ay )£•=!, with ay= / Vr^Vr^j ds. (4.13) 

We also need their scaled versions. Let Dm and be the diagonals of M and A, 
respectively. The scaled matrices are denoted by 

M S = D-/MD~/, A S = D7AD7. (4.14) 

From a simple scaling argument it follows that the spectral condition number of M s 
is bounded uniformly in h and in the shape (ir)rcgularity of the surface triangulation. 
For completeness we include a proof. 
Theorem 4.3. The following holds: 



< 7^ — T < 4 for all v G R JV , v ^ 0. 



v/2 + 2 (DjfV.v 

Proof. The set of all vertices in Fh is denoted by V = { & | 1 < i < N}. Let 

v G l w and vy t G Vu be related by Vh = "i^ii i- e -: v i = v h((,i)- Consider a 

triangle T G J 7 /,, and let its three vertices be denoted by £1,^2, £3- Using quadrature 
we obtain 

Vh{sf ds = ^y(\{vi + v 2 ) 2 + \{V2 + v 3 ) 2 + ^{v 3 + Vl ) 2 ) 

m 

= — (v 1 + v 2 + v 3 + V1V2 + V2V3 + V3V1). 

Hence, J T Vh(s) 2 ds < W~ Y^i=i v 1 holds. From a sign argument it follows that at least 
one of the three terms v\Vi, V2V3 or V3V1 must be positive. Without loss of generality 
we can assume V1V2 > 0. Using \v2v3 + VsVi\ < ^75(^1 + v \ + v i) we get 

v h (s) 2 d8>¥± (v 2 + v 2 + v 2 - ±={v{ + v 2 + v 2 )) = — j3- - (v 2 + v 2 + v 2 ) . 

Note that (Mv,v) = J r Vh(s) 2 ds = ^2 Te:Fh J T fft (s) 2 (is, and thus we obtain, with 
V(T) the set of the three vertices of T, 

£ |T| E ^(0 2 <(Mv,v)< 4l L ^ | T | ^ ^gja. (4.15) 

v + Ter h £ev(T) Te^ h feV(T) 

We observe that 

4 E i t i E ^(0 2 = ^EI su ppWK 2 ( 4 - 16 ) 

Te^ h ?eV(T) i=l 
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holds. From the definition of Dj\/ it follows that 

N „ N 



D M v, v) = W $ ds v\ = £ «? f <?< 

n \T\ 1 N 



f/.S 



(4.17) 



Combination of the results in (|4.15[) . (|4. and (|4.17|) completes the proof. □ 

4.3. Conditioning of the stiffness matrix. We finally address the issue of 
conditioning of the diagonally scaled stiffness matrix A s , cf. (|4.14j) . This matrix has 
a one dimensional kernel due to the constant nodal mode. Thus, we consider the 
effective condition number cond(A s ) = A max (A s )/A2(A s ), where A2 is the minimal 
nonzero eigenvalue. We shall argue below that the condition number of A s can not 
be bounded in general by a constant dependent exclusively on Th, but not on I\. 
Indeed, assume a smooth closed surface T, with |T| = 1, and a smooth function u 
defined on T, such that ||Vru|| i 2( r ) = ||u|| ff 2( r ) = 1. Let Th be the zero level of the 
piecewise linear Lagrange interpolant of the signed distance function to T. Denote 
Uh = IhU e , as in Theorem 14.21 and v — (v%, . . . , v/v) T is the corresponding vector of 
nodal values. From the result in (|4.11l) we obtain 

(Av,v) = ||Vr h « h |Ua(r h ) = l + 0(h). (4.18) 

On the other hand, if there is a node £ in the volume triangulation Th such that 
dist(£, Th) < £ <C 1, then there can appear a triangle in Th with a minimal angle of 
0(e). This implies that there is a diagonal element in A of order 0(e^ 1 ). Without 
lost of generality we may assume An — 0(e~ 1 ) and v\ = 1. Thus we get 

(Dav.v) > A u vl ^Oie- 1 ). (4.19) 

Comparing (|4.18[) and (|4.19[) we conclude that cond(A s ) > 0(e _1 ), with e -> 0. 
Results of numerical experiments in the next section demonstrate that the blow up 
of cond(A s ) can be seen in some cases. 

One might also be interested in a more general dependence of the eigenvalues 
of A s on the distribution of tetrahedral nodes in Td in a neighborhood of IV To a 
certain extend this question is addressed in |S]. 

5. Numerical experiment. In this section we present a few results of numerical 
experiments which illustrate the interpolation estimates from Theorem 14.21 and the 
conditioning of mass and stiffness matrices. Assume the surface T, which is the unit 
sphere r = { x £ M 3 \ \\x\\ = 1 }, is embedded in the bulk domain O = [-2, 2] 3 . The 
signed distance function to T is denoted by d. We construct a hierarchy of uniform 
tetrahedral triangulations {Th} for fi, with h G {1/2, 1/4, 1/8, 1/16, 1/32}. Let d h be 
the piecewise nodal Lagrangian interpolant of d. The triangulated surface is given by 

T h = |J T = {xe n I d h (x) = 0}. 



The corresponding finite clement space Vh consists of all piecewise affinc functions with 
respect to T h , as defined in ([43]) . For h G {1/2,1/4,1/8,1/16,1/32}, the resulting 
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Fig. 5.1. Left: Interpolation error as a function of # d.o.f.; Right: The condition number of 
the mass matrix as a function of of # d.o.f. 

dimensions of Vh arc N = 164,812,3500,14264,57632, respectively. In agreement 
with the 2D nature of T^, we have N ~ hr 2 . 

To illustrate the result of Theorem 14.21 we present the interpolation errors \\u e — 
Ihi^Wtfij-H) an d \u e — IhU e \i t r h for the smooth function 

u(x) = —x\X2 arctan(2a;3) 

defined on the unit sphere, with x = (x\, x%, x^) T . The dependence of the interpola- 
tion errors on the number of degrees of freedom TV is shown in Figure [5TT1 (left). We 
observe the optimal error reduction behavior, consistent with the estimates in (|4.10[) . 

(EU). 

Further, for the same sequence of meshes we compute the spectral condition 
numbers of the mass matrix M and the diagonally scaled mass matrix M s . The 
dependence of the condition numbers on the number of degrees of freedom N is 
illustrated in Figure ISTTI (right). As was proved in Theorem 14. 3[ the scaled mass 
matrix has a uniformly bounded condition number. 

We discussed in section [4731 that concerning the effective condition number of the 
scaled stiffness matrix the situation is more delicate. To illustrate this, we performed 
an experiment in which the intersection between a fixed outer triangulation and the 
surface is varied. Let T be the boundary of the unit sphere with the center located 
in (0, 0, z c ). The discrete surface I\ is defined as described above, induced by the 
uniform outer triangulation. We choose a fixed outer triangulation with h = 1/16. 
We now consider different values for z c , thus "moving the surface through the outer 
triangulation" . The z c values are given in the first column of Table 15.11 Note that for 
the largest shift z c = 0.03 we have z c « 0.5h. 

For the different surface triangulations we computed the interpolation errors as 
described above. It turns out that for the values z c ^ the error behavior is essentially 
the same as that for z c — (illustrated in Fig. I5.ip . 

In the second to fourth columns of Table 15.11 two geometry related quantities 
are given. The second column shows the value of the maximum angle occuring in 
the surface triangulation. Consistent to the theory, cf. Theorem 13. 2\ the maximum 
angle is bounded away from 180°. Small angles, however, can occur. In the third 
and fourth column we show the value of the minimum angle and the number of 
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triangles in the surface triangulation with the smallest angle smaller than 1°. As 
expected, both the minimal angle and this number of small angles strongly varies 
depending on z c . For z c = the smallest angle in the surface triangulation has 
value m in = 1.85°. Extremely small angles can occur, e.g. for z c = 0.00005, we 
have m in = 8.54c-7°. The dimension and the effective condition number of the 
scaled stiffness matrix A s are given in the fifth and sixth column of Table 15.11 The 
values of the condition number show a strong dependence on the sphere location 
(value of z c ). These large condition numbers indicate that linear systems with these 
matrices may be hard to solve using an iterative method. To investigate this further, 
we used the standard PCG MATLAB solver with ILU(0) preconditioner. For given 
v, we computed b = A s v and applied the MATLAB PCG iterative solver with 
a relative residual tolerance of 10 -8 . The resulting iteration numbers are given in 
column 7 of Tabic 15.11 These iteration counts are "high" compared to the ones 
that are generally needed for standard discretization of diffusion problems. To make 
this more quantitative, we constructed a reference matrix A rof as follows: A rcf = 
blocktridiag(-B T , D, -B), with D = tridiag(-l, 6,-1), B = tridiag(0, 1, 1). In most 
rows the matrix A has 7 nonzero entries, which is approximately the same as the 
average number of nonzero entries per row in the matrix A. s used in the experiment. In 
A ref we use 120 blockrows and blockcolums and the matrices D and B have dimension 
120. Then the matrix A has dimension 14400, which is comparable to the dimension 
of A s used in the experiment, cf. Table HTT1 The same iterative solver with the same 
stopping criterion applied to a linear system with A resulted in 42 PCG iterations, 
which is much lower than the iteration numbers listed in Tabic 15.11 
In view of these observations, and the fact that solving a PDE on a surface (in 
3D) is a iwo-dimensional problem, it is better to use a direct solver. We performed 
experiments with the MATLAB sparse direct solver A s \ b. We measured computing 
time by the MATLAB function CPUTIME. For the system with the reference matrix 
A ref we obtained (on our machine) CPUTIME= 1.38. For the matrix A s we obtained 
CPU time measurements given in the last column of Table 15.11 These show that for 
the direct MATLAB solver the matrices A s are not (much) more difficult to deal 
with than the reference matrix A rcf . Variations in CPU times are probably caused by 
slightly different fill-in properties of matrices for different grids. The one dimensional 
kernel of the matrix A s did not cause difficulties for the solver. We checked the 
accuracy of the computed solution (in the energy norm) and this was satisfactory. 
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0m ax 


0min 


#T: 
<*min<l° 


dim(A s ) 


cond(A s ) 


# PCG 


CPUTIME 


0.03 


147.4° 


0.050° 


420 


14406 


1.82e+4 


245 


3.64 


0.02 


145.3° 


0.027° 


292 


14376 


2.20e+4 


282 


3.52 


0.008 


145.4° 


0.014° 


270 


14368 


3.44e+4 


331 


3.61 


0.002 


144.3° 


0.002° 


126 


14300 


1.94e+5 


285 


2.33 


0.0005 


141.0° 


1.22e-4° 


20 


14288 


3.07e+6 


259 


1.93 


0.00025 


140.4° 


3.05e-5° 


20 


14288 


1.23e+7 


191 


2.22 


0.00005 


139.9° 


8.54e-7° 


24 


14288 


3.06e+8 


202 


1.43 





139.8° 


1.85° 





14264 


9.14e+3 


142 


2.85 



Table 5.1 

Angles in the surface triangulation, dimension of A 3 , cond(A s ), iteration count for PCG and 
timing for direct solver. 



6. Conclusions. The main new result of this paper is a geometric property of 
the piecewise planar surface which is the zero level of a continuous piecewise affinc 
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level set function. If this pieccwisc planar surface is consistent with an outer tetra- 
hcdral triangulation that satisfies the minimum angle condition, then after a suitable 
subdivision of the quadrilaterals into two triangles the resulting surface triangulation 
satisfies a maximum angle condition. This maximum angle property of the surface 
triangulation is used to derive optimal error bounds for the nodal interpolation oper- 
ator in the finite element space of continuous pieccwisc linear functions on the surface 
triangulation. This implies that the discretization of a surface diffusion PDE in this 
finite element space results in optimal discretization error bounds. We study the 
conditioning of the scaled mass and stiffness matrices corresponding to this finite 
element space. The condition number of the scaled mass matrix is shown to be uni- 
formly bounded. The scaled stiffness matrix can have a very large effective condition 
number. Results of a numerical experiment indicate that for solving systems with the 
scaled stiffness matrix it is better to use a sparse direct solver rather than an iterative 
solver. A topic that we plan to investigate further is whether some grid smoothing 
(elimination of extremely small angles) can be developed such that the optimal ap- 
proximation property still holds and the conditioning of the scaled stiffness matrix is 
improved. 
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